function [tOut, uOut, udtOut, udtdtOut] = h2hydro(params)

  Nx = params.Nx;
  Ny = params.Ny;


  % ****************************
  % initial conditions
  % ****************************
  T  = params.TempOffset .* ones(Nx, Ny);
  ux = zeros(Nx, Ny) + 0.0001;
  uy = zeros(Nx, Ny) + 0.0001;
  P1 = zeros(Nx, Ny);
  P2 = zeros(Nx, Ny);

  % ********************************************************
  %     V O R T I C E S
  % ********************************************************

  % initial conditions 
  if params.hydroVortices
    vp = params.hydroVorticesParams;
    for i =1:size(vp,1)
      v = vp(i,:);
      vamp = v(4);
      vwidth = v(3);
      vx = v(1);
      vy = v(2);

      cx = vx.*params.Nx;
      cy = vy.*params.Ny;

      gauss_grid_width = vwidth .* params.Nx;
      sig2 = 2.*(gauss_grid_width).^2;

      for x = 1:params.Nx
        for y = 1:params.Nx
          for ddx = -1:1
            for ddy = -1:1
              %T(x, y) = T(x,y) + params.hydroPeakAmp .* exp(-((x-cx)^2 + (y-cy)^2)/sig2);
              ux(x,y) = ux(x,y) -vamp.*(y+ddy*params.Ny-cy).*exp(-((x+ddx*params.Nx-cx)^2 + (y+ddy*params.Ny-cy)^2)/sig2);
              uy(x,y) = uy(x,y) +vamp.*(x+ddx*params.Nx-cx).*exp(-((x+ddx*params.Nx-cx)^2 + (y+ddy*params.Ny-cy)^2)/sig2);
            end
          end
        end
      end  % end of fors
    end
  end

  % ********************************************************
  %     E X P    P E A K
  % ********************************************************

  if params.hydroPeak
    cx = 0.5 * params.Nx + 0.5;
    cy = 0.5 * params.Ny + 0.5;
    cx2 = 0.5 * params.Nx + 0.5;
    cy2 = 0.6 * params.Ny + 0.5;
    gauss_grid_width = params.hydroPeakRelWidth .* params.Nx;

    sig2 = 2.*(gauss_grid_width).^2;
    for x = 1:params.Nx
      for y = 1:params.Nx
        T(x, y) = T(x,y) + params.hydroPeakAmp .* exp(-((x-cx)^2 + (y-cy)^2)/sig2);
      end
    end
  end
  %ux = 0.99999999999999;

  [gx, gy] = ndgrid(params.Cx/params.Lx, params.Cy/params.Ly);
  gx = 2.*pi.*gx;
  gy = 2.*pi.*gy;


  % ********************************************************
  %     R A N D O M
  % ********************************************************

    rng(123);
  if params.hydroRandGradOn
    % smooth random gradients
    rng(123);
    N = 1000;
    %random_factor = 1/50; %correct for random interference
    for s=1:N
      phasex = rand(1)*2*pi;
      phasey = rand(1)*2*pi;
      ampx = rand(1);
      ampy = rand(1);
      freqx  = round(rand(1)*params.hydroRandGradFreq);
      freqy  = round(rand(1)*params.hydroRandGradFreq);
      ux = ux + 1/N.*ampx.*sin(freqx.*gx + phasex).*sin(freqy.*gy + phasey);

      phasex = rand(1)*2*pi;
      phasey = rand(1)*2*pi;
      ampx = rand(1)*params.hydroPureWaveAmp;
      freqx  = round(rand(1)*params.hydroRandGradFreq);
      freqy  = round(rand(1)*params.hydroRandGradFreq);
      uy = uy + 1/N.*ampx.*sin(freqx.*gx + phasex).*sin(freqy.*gy + phasey); %.*sin(freqy*gy + phasey);
    end
    maxx = max(max(abs(ux)));
    maxy = max(max(abs(ux)));
    disp(sprintf('max ux: %f  max uy: %f', maxx, maxy));

    if true
      normalization = 0.1;
      ux = ux / max([maxx maxy]) .* normalization;
      uy = uy / max([maxx maxy]) .* normalization;
    end
  end

    
  if true
    keyboard()
    To = 1;   % Temp = (3/4pi)T = 3/4pi
    Nw = 10;  % Nw full oscillations of uy in x
    Nk = 1;   % Number of fourier modes contributing to noise
    du = 1/5;
    dux = 0*X;
    duy = 0*X;
    for kx = 1:Nk
        for ky = 1:Nk
            dux = dux + rand(1) * cos( kx.* X + ky.* Y + 2*pi*rand(1));
            duy = duy + rand(1) * cos( kx.* X + ky.* Y + 2*pi*rand(1));
        end
    end
    dux = dux / (max(dux(:))/du);
    duy = duy / (max(duy(:))/du);


    ux =                      dux;
    uy = cos((Nw*pi/Lx)*Xp) + duy;

    % turbulent flow
    % ux = ux + 0.1 * sin(2*gy);%.*(1+0.2.*cos(2*gx));
    % uy = uy + 0.005.*sin(2*gx+rand(1)*cos(3*gy))
    % uy = uy + 0.003.*sin(1*gx+rand(1)*cos(2*gy))
    % uy = uy + 0.0002.*sin(5*gx+rand(1)*cos(4*gx))
    % amp = 0.1;
    % P1  = P1 + amp*1*sin(4*gy+rand(1));
    % P2  = P2 + amp*1*sin(3*gx+rand(1));
    ux = 
  end


  if params.hydroPureWave
    %ux = ux + params.hydroPureWaveAmp.*sin(params.hydroPureWaveFreq.*gy);
  end

  if false
    uxmax = max(max(abs(ux)));
    uymax = max(max(abs(uy)));
    umax = max([uxmax uymax]);
    ux = ux./umax * 0.5;
    uy = uy./umax * 0.5;
  end

  if false
    % copied from allan's file
    x = -pi + 2*pi*(0:Nx-1)'/Nx;
    y = -pi + 2*pi*(0:Ny-1)'/Ny;
    [X,Y]=ndgrid(x,y);

    % T = 0.*T + 1.0/50;
    ux = 0*X + 0.5 * sin(2*Y);
    uy = 0*X + 0.0001*sin(3*X+rand(1)*cos(7*Y));
    P1  = 0.000001*sin(3*X+rand(1));
    P2  = 0.000001*sin(2*Y+rand(1));



      cx = 0.5 .* params.Nx;
      cy = 0.3 .* params.Ny;
      sig2 = (0.05 .* params.Nx).^2 * 2;
      for x = 1:params.Nx
        for y = 1:params.Nx
          for ddx = -1:1
            for ddy = -1:1
              uy(x, y) = uy(x,y) + 0.05.* exp(-((x-cx)^2 + (y-cy)^2)/sig2);
              %ux(x,y) = ux(x,y) + rand*0.0001;;
              %uy(x,y) = uy(x,y) + rand*0.0001;;
            end
          end
        end
      end  % end of fors
  end

  % noise

  % calculate missing initial conditions
  %[P1 P2] = h2icgen(T, ux, uy, params);


  % ****************************
  %  rescale temperature
  % ****************************

  % T = T .* (4 * pi / 3);

  % repackage initial conditions into a vector
  Nt = Nx.*Ny;
  uIn = zeros(Nt.*5, 1);
  uIn(0*Nt+1 : 1*Nt) = T(:);
  uIn(1*Nt+1 : 2*Nt) = ux(:);
  uIn(2*Nt+1 : 3*Nt) = uy(:);
  uIn(3*Nt+1 : 4*Nt) = P1(:);
  uIn(4*Nt+1 : 5*Nt) = P2(:);


  % ****************************
  % run the ode solver
  % ****************************

  time_steps = (0:params.Tsteps-1)/params.Tsteps*params.Tfinal;

  tic;
  %uOut   = zeros(params.Tsteps, Nt.*5);
  %tOut   = zeros(params.Tsteps, Nt.*5);

    opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-8, 'OutputFcn', @(t,u,flag) h2out(t, u, flag, params));

  if params.hydroOrder == 2
    [tOut, uOut] = ode45(@(t,u) h2eom(t,u,params), time_steps, uIn, opts);
  end
  if params.hydroOrder == -2
    [tOut, uOut] = ode45(@(t,u) h2eomnr(t,u,params), time_steps, uIn, opts);
  end
  if params.hydroOrder == -3
    [tOut, uOut] = ode45(@(t,u) h2eomnr2(t,u,params), time_steps, uIn, opts);
  end
  if params.hydroOrder == 0
    [tOut, uOut] = ode45(@(t,u) h0eom(t,u,params), time_steps, uIn, opts);
  end

  elapsed = toc;
  disp(sprintf('computed hydro simulation in %fsec', elapsed));


  % ****************************
  % compute 1st time derivatives
  % ****************************

  tic;

  udtOut   = zeros(params.Tsteps, Nt.*5);

  % first derivative from eom
  for it=1:params.Tsteps
    t = tOut(it);
    u = uOut(it, :);
    
    udtOut(it,:) = h2eom(t, u, params);
  end

  % ****************************
  %  get rid of Pi and B
  % ****************************

  uOut  (:,3*Nt+1 : 5*Nt) = [];
  udtOut(:,3*Nt+1 : 5*Nt) = [];

  % ****************************
  % compute 2nd time derivatives
  % ****************************

  derorder = params.timeOrder;
  % second der from finite differences, O(8)
  derker = d1ker(params.timeOrder);
  %derker = [1.0/12; -2.0/3; 0; 2.0/3; -1.0/12]; % o4
  udtdtOut = conv2(udtOut, derker, 'same')./params.Dt;

  % delete derivative padding
  deli = [1:(derorder/2) (length(tOut)-derorder/2+1):length(tOut)];
  tOut    (deli, :) = [];
  uOut    (deli, :) = [];
  udtOut  (deli, :) = [];
  udtdtOut(deli, :) = [];
  elapsed = toc;
  disp(sprintf('computed higher hydro time derivatives in %fsec', elapsed));

  params.Tbegin = params.Tbegin + derorder./2.*params.Dt;
  params.Tfinal = params.Tfinal - derorder./2.*params.Dt;
  params.Tsteps = params.Tsteps - derorder;


  % ****************************
  % compute ut
  % ****************************

  ux = uOut(:, 1*Nt+1 : 2*Nt);
  uy = uOut(:, 2*Nt+1 : 3*Nt);

  uxdt = udtOut(:, 1*Nt+1 : 2*Nt);
  uydt = udtOut(:, 2*Nt+1 : 3*Nt);

  uxdtdt = udtdtOut(:, 1*Nt+1 : 2*Nt);
  uydtdt = udtdtOut(:, 2*Nt+1 : 3*Nt);

  ut = sqrt(1+ux.^2+uy.^2);
  utdt = (uxdt.*ux + uydt.*uy)./ut;
  utdtdt = (1+ux.^2+uy.^2).^(-3/2).*(uxdt.^2.*(1+uy.^2)+(-2).*ux.*uxdt.*uy.*uydt+(1+ux.^2).*uydt.^2+(1+ux.^2+uy.^2).*(ux.*uxdtdt+uy.*uydtdt));

  uOut    (:, 3*Nt+1 : 4*Nt) = ut;
  udtOut  (:, 3*Nt+1 : 4*Nt) = utdt;
  udtdtOut(:, 3*Nt+1 : 4*Nt) = utdtdt;
  
  % can safely exit here if you want to ignore all the stress energy computations
  return
  if (false)
  for ti = 1:size(uOut,1)
    T  = reshape(uOut(ti, 0*Nt+1 : 1*Nt), Nx, Ny);
    ux = reshape(uOut(ti, 1*Nt+1 : 2*Nt), Nx, Ny);
    uy = reshape(uOut(ti, 2*Nt+1 : 3*Nt), Nx, Ny);
    ut = reshape(uOut(ti, 3*Nt+1 : 4*Nt), Nx, Ny);


     Tdt = reshape(udtOut(ti, 0*Nt+1 : 1*Nt), Nx, Ny);
    uxdt = reshape(udtOut(ti, 1*Nt+1 : 2*Nt), Nx, Ny);
    uydt = reshape(udtOut(ti, 2*Nt+1 : 3*Nt), Nx, Ny);
    utdt = reshape(udtOut(ti, 3*Nt+1 : 4*Nt), Nx, Ny);

    uxdtdt = reshape(udtdtOut(ti, 1*Nt+1 : 2*Nt), Nx, Ny);
    uydtdt = reshape(udtdtOut(ti, 2*Nt+1 : 3*Nt), Nx, Ny);
    utdtdt = reshape(udtdtOut(ti, 3*Nt+1 : 4*Nt), Nx, Ny);

    Tdx = dx(T);
    Tdy = dy(T);

    utdx = dx(ut);
    utdy = dy(ut);

    uxdx = dx(ux);
    uxdy = dy(ux);

    uydx = dx(uy);
    uydy = dy(uy);

    utdtdx  = dx(utdt);
    utdtdy  = dy(utdt);
    utdxdx  = dx(utdx);
    utdxdy  = dy(utdx);
    utdydy  = dy(utdy);

    uxdtdx  = dx(uxdt);
    uxdtdy  = dy(uxdt);
    uxdxdx  = dx(uxdx);
    uxdxdy  = dy(uxdx);
    uxdydy  = dy(uxdy);

    uydtdx  = dx(uydt);
    uydtdy  = dy(uydt);
    uydxdx  = dx(uydx);
    uydxdy  = dy(uydx);
    uydydy  = dy(uydy);
    
    % d_mu T^mu 0
    eomo1 = zeros(3, params.Nx, params.Ny);
    eomo1(1,:,:) = (1/2).*T.*(3.*T.^2.*(utdx.*ux+utdy.*uy+ut.*(2.*utdt+uxdx+uydy))+2.*(Tdt.*utdt+(-2).*Tdt.*ut.^2.*utdt+Tdt.*ut.^4.*utdt+(-1).*Tdx.*utdx+Tdx.*ut.^2.*utdx+Tdx.*uxdt+(-1).*Tdx.*ut.^2.*uxdt+(-1).*Tdt.*uxdx+Tdt.*ut.^2.*uxdx+(-1).*Tdx.*ut.*ux.^3.*uxdx+(-1).*Tdt.*ut.*utdy.*uy+Tdt.*ut.^3.*utdy.*uy+(-1).*Tdx.*ut.*uxdy.*uy+Tdt.*ut.*uy.*uydt+(-1).*Tdt.*ut.^3.*uy.*uydt+(-1).*Tdx.*ut.*uy.*uydx+ux.^2.*(Tdx.*((-1).*utdx+uxdt)+(-1).*Tdt.*uxdx+ut.^2.*(Tdx.*(utdx+(-1).*uxdt)+(-1).*Tdt.*uxdx)+(-1).*Tdx.*ut.*uy.*(uxdy+uydx))+(-1).*Tdt.*uydy+Tdt.*ut.^2.*uydy+(-1).*Tdt.*uy.^2.*uydy+(-1).*Tdt.*ut.^2.*uy.^2.*uydy+(-1).*ux.*(uy.*((-1).*Tdx.*((-1)+ut.^2).*utdy+Tdt.*(1+ut.^2).*uxdy+(-1).*Tdx.*uydt+Tdx.*ut.^2.*uydt+Tdt.*uydx+Tdt.*ut.^2.*uydx)+Tdx.*ut.*uy.^2.*uydy+ut.*((-1).*Tdt.*((-1)+ut.^2).*(utdx+(-1).*uxdt)+Tdx.*((-1).*((-1)+ut.^2).*utdt+uxdx)+(-1).*Tdx.*uydy))+Tdy.*(ut.^3.*utdt.*uy+(-1).*utdx.*ux.*uy+ux.*uxdt.*uy+((-1)+ut.^2).*utdy.*(1+uy.^2)+uydt+uy.^2.*uydt+ut.^2.*(ux.*(utdx+(-1).*uxdt).*uy+(-1).*(1+uy.^2).*uydt)+(-1).*ut.*(ux.^2.*uxdx.*uy+ux.*(1+uy.^2).*(uxdy+uydx)+uy.*(utdt+(-1).*uxdx+(1+uy.^2).*uydy))))+T.*((-3).*Tdt+utdtdt+ut.^4.*utdtdt+(-1).*utdxdx+(-1).*utdydy+(-1).*utdx.*ux.^3.*uxdx+(-2).*utdt.*utdy.*uy+(-2).*utdx.*uxdy.*uy+(-1).*utdydy.*uy.^2+utdt.*uy.*uydt+uxdx.*uy.*uydt+(-1).*utdx.*uy.*uydx+(-1).*uxdt.*uy.*uydx+(-1).*ux.^2.*(utdxdx+uy.*(utdy.*uxdx+utdx.*uxdy+utdx.*uydx))+(-3).*utdy.*uy.*uydy+(-1).*utdy.*uy.^3.*uydy+ut.^3.*(4.*utdt.^2+2.*utdtdx.*ux+utdx.*uxdt+(-1).*uxdt.^2+(-1).*ux.*uxdtdt+utdt.*uxdx+2.*utdtdy.*uy+utdy.*uydt+(-1).*uydt.^2+(-1).*uy.*uydtdt+utdt.*uydy)+(-1).*ux.*(2.*utdt.*utdx+(-1).*utdt.*uxdt+3.*utdx.*uxdx+2.*utdxdy.*uy+uxdy.*uydt+utdy.*(uxdy.*(1+uy.^2)+(2+uy.^2).*uydx)+(-1).*uxdt.*uydy+utdx.*uy.^2.*uydy)+ut.^2.*(9.*Tdt+(-2).*utdtdt+utdxdx+utdydy+6.*utdt.*utdx.*ux+utdxdx.*ux.^2+(-3).*utdt.*ux.*uxdt+(-2).*ux.^2.*uxdtdx+2.*utdx.*ux.*uxdx+(-4).*ux.*uxdt.*uxdx+(-1).*ux.*uxdy.*uydt+uy.^2.*(utdydy+(-2).*uydtdy)+utdy.*ux.*uydx+(-2).*ux.*uydt.*uydx+utdx.*ux.*uydy+(-1).*ux.*uxdt.*uydy+uy.*(2.*utdxdy.*ux+(-2).*ux.*uxdtdy+(utdx+(-2).*uxdt).*uxdy+(-3).*utdt.*uydt+(-1).*uxdx.*uydt+(-2).*ux.*uydtdx+(-1).*uxdt.*uydx+(-4).*uydt.*uydy+utdy.*(6.*utdt+uxdx+2.*uydy)))+(-1).*ut.*(4.*utdt.^2+(-2).*utdx.^2+(-2).*utdy.^2+(-9).*Tdx.*ux+2.*utdtdx.*ux+(-2).*utdx.^2.*ux.^2+3.*utdx.*uxdt+2.*utdx.*ux.^2.*uxdt+(-1).*uxdt.^2+(-1).*ux.*uxdtdt+(-1).*utdt.*uxdx+2.*utdt.*ux.^2.*uxdx+uxdx.^2+3.*ux.^2.*uxdx.^2+ux.*uxdxdx+ux.^3.*uxdxdx+uxdy.^2+ux.*uxdydy+3.*utdy.*uydt+(-1).*uydt.^2+2.*uxdy.*uydx+ux.^2.*uxdy.*uydx+uydx.^2+ux.^2.*uydx.^2+(-1).*utdt.*uydy+(-2).*uxdx.*uydy+ux.^2.*uxdx.*uydy+uydy.^2+uy.^2.*((-2).*utdy.^2+uxdy.^2+ux.*uxdydy+2.*utdy.*uydt+uxdy.*uydx+2.*ux.*uydxdy+2.*utdt.*uydy+uxdx.*uydy+3.*uydy.^2)+uy.^3.*uydydy+uy.*((-9).*Tdy+2.*utdtdy+(-1).*uydtdt+uydxdx+ux.^2.*(2.*uxdxdy+uydxdx)+2.*ux.*(utdy.*((-2).*utdx+uxdt)+utdx.*uydt+utdt.*uydx+uxdx.*uydx+2.*uydx.*uydy+uxdy.*(utdt+2.*uxdx+uydy))+uydydy))));
     eomo1(2,:,:) = (1/2).*T.*(3.*T.^2.*(ut.*uxdt+uxdy.*uy+ux.*(utdt+2.*uxdx+uydy))+2.*(Tdx.*utdt+(-1).*Tdt.*utdx+Tdt.*ut.^3.*utdt.*ux+Tdx.*utdt.*ux.^2+(-1).*Tdt.*utdx.*ux.^2+Tdt.*uxdt+Tdt.*ux.^2.*uxdt+(-1).*Tdx.*uxdx+(-2).*Tdx.*ux.^2.*uxdx+(-1).*Tdx.*ux.^4.*uxdx+(-1).*Tdt.*utdy.*ux.*uy+(-1).*Tdx.*ux.*uxdy.*uy+(-1).*Tdx.*ux.^3.*uxdy.*uy+Tdt.*ux.*uy.*uydt+ut.^2.*((-1).*Tdx.*utdt+Tdt.*utdx+(-1).*Tdt.*uxdt+ux.^2.*(Tdx.*utdt+Tdt.*utdx+(-1).*Tdt.*uxdt)+ux.*uy.*(Tdy.*utdt+Tdt.*utdy+(-1).*Tdt.*uydt))+(-1).*Tdx.*ux.*uy.*uydx+(-1).*Tdx.*ux.^3.*uy.*uydx+Tdx.*uydy+Tdx.*ux.^2.*uydy+Tdx.*uy.^2.*uydy+(-1).*Tdx.*ux.^2.*uy.^2.*uydy+(-1).*Tdy.*(ux.^3.*uxdx.*uy+(1+ux.^2).*uxdy.*(1+uy.^2)+(1+uy.^2).*uydx+ux.^2.*(1+uy.^2).*uydx+ux.*uy.*((-1).*utdt+uxdx+(1+uy.^2).*uydy))+ut.*(ux.^3.*(Tdx.*(utdx+(-1).*uxdt)+(-1).*Tdt.*uxdx)+ux.^2.*uy.*(Tdx.*utdy+Tdy.*(utdx+(-1).*uxdt)+(-1).*Tdt.*uxdy+(-1).*Tdx.*uydt+(-1).*Tdt.*uydx)+(-1).*uy.*(Tdx.*utdy+Tdy.*((-1).*utdx+uxdt)+Tdt.*uxdy+(-1).*Tdx.*uydt+Tdt.*uydx)+ux.*(Tdx.*(utdx+(-1).*uxdt)+Tdy.*(1+uy.^2).*(utdy+(-1).*uydt)+(-1).*Tdt.*(utdt+uxdx+((-1)+uy.^2).*uydy))))+(-1).*T.*((-3).*Tdx+(-1).*ut.^3.*utdt.*uxdt+(-1).*uxdtdt+uxdxdx+ux.^4.*uxdxdx+uxdydy+2.*utdy.*uxdt.*uy+2.*uxdx.*uxdy.*uy+uxdydy.*uy.^2+ut.^2.*(uxdtdt+(-1).*uy.*(utdt.*uxdy+uxdt.*(utdy+(-1).*uydt)))+(-1).*utdx.*uy.*uydt+(-1).*uxdt.*uy.*uydt+utdt.*uy.*uydx+uxdx.*uy.*uydx+3.*uxdy.*uy.*uydy+uxdy.*uy.^3.*uydy+ux.^3.*((-1).*utdx.^2+(-1).*ut.*utdxdx+utdx.*uxdt+2.*ut.*uxdtdx+utdt.*uxdx+4.*uxdx.^2+2.*uxdxdy.*uy+uxdy.*uydx+uydx.^2+uy.*uydxdx+uxdx.*uydy)+ut.*(3.*utdt.*uxdt+(-1).*utdx.*uxdx+2.*uxdt.*uxdx+2.*uxdtdy.*uy+2.*uxdy.*uydt+uxdy.*uy.^2.*uydt+(-1).*utdy.*(uxdy.*(1+uy.^2)+(-1).*uydx)+((-1).*utdx+uxdt.*uy.^2).*uydy)+ux.^2.*((-9).*Tdx+(-2).*ut.^2.*utdtdx+(-4).*ut.*utdt.*utdx+2.*ut.*utdt.*uxdt+(-1).*uxdtdt+ut.^2.*uxdtdt+(-3).*ut.*utdx.*uxdx+6.*ut.*uxdt.*uxdx+2.*uxdxdx+uxdydy+ut.*uxdy.*uydt+(-1).*ut.*utdy.*uydx+2.*ut.*uydt.*uydx+uy.^2.*(uxdydy+2.*uydxdy)+(-1).*ut.*utdx.*uydy+ut.*uxdt.*uydy+uy.*((-2).*ut.*utdxdy+utdy.*((-2).*utdx+uxdt)+2.*ut.*uxdtdy+utdx.*uydt+2.*ut.*uydtdx+utdt.*uydx+3.*uxdx.*uydx+4.*uydx.*uydy+uxdy.*(utdt+6.*uxdx+2.*uydy)))+ux.*((-9).*Tdt.*ut+utdt.^2+(-3).*ut.^2.*utdt.^2+ut.*utdtdt+(-1).*ut.^3.*utdtdt+(-1).*utdx.^2+(-1).*ut.*utdxdx+(-1).*utdy.^2+(-1).*ut.*utdydy+3.*utdx.*uxdt+(-2).*ut.^2.*utdx.*uxdt+(-2).*uxdt.^2+2.*ut.^2.*uxdt.^2+2.*ut.*uxdtdx+(-1).*utdt.*uxdx+(-2).*ut.^2.*utdt.*uxdx+4.*uxdx.^2+2.*uxdy.^2+2.*utdy.*uydt+(-1).*ut.^2.*utdy.*uydt+(-1).*uydt.^2+ut.^2.*uydt.^2+3.*uxdy.*uydx+uydx.^2+(-2).*utdt.*uydy+(-1).*ut.^2.*utdt.*uydy+(-1).*uxdx.*uydy+uydy.^2+uy.^2.*((-1).*utdy.^2+(-1).*ut.*utdydy+2.*uxdy.^2+utdy.*uydt+2.*ut.*uydtdy+2.*uxdy.*uydx+utdt.*uydy+2.*uxdx.*uydy+3.*uydy.^2)+uy.^3.*uydydy+uy.*((-9).*Tdy+2.*uxdxdy+(-1).*uydtdt+ut.^2.*((-2).*utdtdy+uydtdt)+uydxdx+(-2).*ut.*((utdx+(-2).*uxdt).*uxdy+(-1).*utdt.*uydt+(-1).*uxdx.*uydt+(-1).*uxdt.*uydx+(-2).*uydt.*uydy+utdy.*(2.*utdt+uxdx+uydy))+uydydy))));
     eomo1(3,:,:) = (1/2).*T.*(3.*T.^2.*(ut.*uydt+ux.*uydx+uy.*(utdt+uxdx+2.*uydy))+2.*(Tdt.*ut.^3.*utdt.*uy+(-1).*Tdx.*ux.^3.*uxdx.*uy+ut.^2.*((-1).*Tdy.*utdt+ux.*(Tdx.*utdt+Tdt.*utdx+(-1).*Tdt.*uxdt).*uy+Tdy.*utdt.*uy.^2+Tdt.*utdy.*(1+uy.^2)+(-1).*Tdt.*uydt+(-1).*Tdt.*uy.^2.*uydt)+(-1).*ux.^2.*(Tdy.*uxdx.*((-1)+uy.^2)+Tdx.*uxdy.*(1+uy.^2)+Tdx.*(1+uy.^2).*uydx)+(-1).*ux.*uy.*((-1).*Tdx.*utdt+Tdt.*utdx+(-1).*Tdt.*uxdt+Tdx.*uxdx+Tdy.*(1+uy.^2).*(uxdy+uydx)+Tdx.*(1+uy.^2).*uydy)+(-1).*(1+uy.^2).*(Tdt.*utdy+Tdx.*uxdy+(-1).*Tdt.*uydt+Tdx.*uydx+Tdy.*((-1).*utdt+(-1).*uxdx+(1+uy.^2).*uydy))+ut.*(ux.*uy.^2.*(Tdx.*utdy+Tdy.*(utdx+(-1).*uxdt)+(-1).*Tdt.*uxdy+(-1).*Tdx.*uydt+(-1).*Tdt.*uydx)+ux.*(Tdx.*utdy+Tdy.*((-1).*utdx+uxdt)+(-1).*Tdt.*uxdy+(-1).*Tdx.*uydt+(-1).*Tdt.*uydx)+uy.^3.*(Tdy.*(utdy+(-1).*uydt)+(-1).*Tdt.*uydy)+uy.*(Tdx.*(1+ux.^2).*(utdx+(-1).*uxdt)+Tdy.*(utdy+(-1).*uydt)+(-1).*Tdt.*(utdt+((-1)+ux.^2).*uxdx+uydy))))+(-1).*T.*((-1).*utdy.*ux.*uxdt+(-1).*ut.*utdy.*uxdx+ut.*utdx.*uxdy+utdt.*ux.*uxdy+(-3).*Tdy.*(1+3.*uy.^2)+3.*ut.*utdt.*uydt+(-1).*ut.^3.*utdt.*uydt+2.*utdx.*ux.*uydt+(-1).*ut.^2.*utdx.*ux.*uydt+(-1).*ux.*uxdt.*uydt+ut.^2.*ux.*uxdt.*uydt+ut.*ux.^2.*uxdx.*uydt+(-1).*uydtdt+ut.^2.*uydtdt+2.*ut.*ux.*uydtdx+(-1).*ut.*utdx.*uydx+(-1).*ut.^2.*utdt.*ux.*uydx+(-1).*ut.*utdx.*ux.^2.*uydx+2.*ut.*uxdt.*uydx+ut.*ux.^2.*uxdt.*uydx+3.*ux.*uxdx.*uydx+ux.^3.*uxdx.*uydx+uydxdx+ux.^2.*uydxdx+(-1).*ut.*utdy.*uydy+ux.*uxdy.*uydy+2.*ut.*uydt.*uydy+2.*ux.*uydx.*uydy+uy.^3.*((-1).*utdy.^2+(-1).*ut.*utdydy+uxdy.^2+ux.*uxdydy+utdy.*uydt+2.*ut.*uydtdy+uxdy.*uydx+2.*ux.*uydxdy+utdt.*uydy+uxdx.*uydy+4.*uydy.^2)+uy.*((-9).*Tdt.*ut+utdt.^2+(-3).*ut.^2.*utdt.^2+ut.*utdtdt+(-1).*ut.^3.*utdtdt+(-1).*utdx.^2+(-1).*ut.*utdxdx+(-1).*utdy.^2+(-1).*ut.*utdydy+(-9).*Tdx.*ux+(-2).*ut.^2.*utdtdx.*ux+(-4).*ut.*utdt.*utdx.*ux+(-1).*utdx.^2.*ux.^2+(-1).*ut.*utdxdx.*ux.^2+2.*utdx.*uxdt+(-1).*ut.^2.*utdx.*uxdt+2.*ut.*utdt.*ux.*uxdt+utdx.*ux.^2.*uxdt+(-1).*uxdt.^2+ut.^2.*uxdt.^2+(-1).*ux.*uxdtdt+ut.^2.*ux.*uxdtdt+2.*ut.*ux.^2.*uxdtdx+(-2).*utdt.*uxdx+(-1).*ut.^2.*utdt.*uxdx+(-2).*ut.*utdx.*ux.*uxdx+utdt.*ux.^2.*uxdx+4.*ut.*ux.*uxdt.*uxdx+uxdx.^2+3.*ux.^2.*uxdx.^2+ux.*uxdxdx+ux.^3.*uxdxdx+uxdy.^2+ux.*uxdydy+(-2).*uydt.^2+2.*ut.^2.*uydt.^2+2.*ut.*uydtdy+4.*ut.*ux.*uydt.*uydx+2.*uydx.^2+2.*ux.^2.*uydx.^2+utdy.*(3.*uydt+(-2).*ut.^2.*uydt+(-2).*ut.*ux.*uydx)+uxdy.*(2.*ut.*ux.*uydt+(3+2.*ux.^2).*uydx)+2.*ux.*uydxdy+(-1).*utdt.*uydy+(-2).*ut.^2.*utdt.*uydy+(-2).*ut.*utdx.*ux.*uydy+2.*ut.*ux.*uxdt.*uydy+(-1).*uxdx.*uydy+2.*ux.^2.*uxdx.*uydy+4.*uydy.^2)+uydydy+uy.^4.*uydydy+uy.^2.*((-1).*uydtdt+ut.^2.*((-2).*utdtdy+uydtdt)+uydxdx+ux.^2.*(2.*uxdxdy+uydxdx)+ut.*((-2).*utdxdy.*ux+2.*ux.*uxdtdy+(-1).*(utdx+(-2).*uxdt).*uxdy+2.*utdt.*uydt+uxdx.*uydt+2.*ux.*uydtdx+uxdt.*uydx+6.*uydt.*uydy+(-1).*utdy.*(4.*utdt+uxdx+3.*uydy))+ux.*(utdy.*((-2).*utdx+uxdt)+utdx.*uydt+utdt.*uydx+2.*uxdx.*uydx+6.*uydx.*uydy+uxdy.*(utdt+4.*uxdx+3.*uydy))+2.*uydydy)));
     eomo0(1,:,:) = (3/2).*T.^2.*((-1).*Tdt+3.*Tdt.*ut.^2+T.*(utdx.*ux+utdy.*uy)+ut.*(3.*Tdx.*ux+3.*Tdy.*uy+T.*(2.*utdt+uxdx+uydy)));
     eomo0(2,:,:) = (3/2).*T.^2.*(Tdx+3.*Tdx.*ux.^2+T.*(ut.*uxdt+uxdy.*uy)+ux.*(3.*Tdt.*ut+3.*Tdy.*uy+T.*(utdt+2.*uxdx+uydy)));
     eomo0(3,:,:) = (3/2).*T.^2.*(3.*(Tdt.*ut+Tdx.*ux).*uy+Tdy.*(1+3.*uy.^2)+T.*(ut.*uydt+ux.*uydx+uy.*(utdt+uxdx+2.*uydy)));
     
     for  i=1:3
       erravg(ti, 1,i) = sum(sum(abs(eomo0(i,:,:)))).*params.Dx.*params.Dy/params.Lx/params.Ly;
       errmax(ti, 1,i) = max(max(abs(eomo0(i,:,:))));
       erravg(ti, 2,i) = sum(sum(abs(eomo1(i,:,:)))).*params.Dx.*params.Dy/params.Lx/params.Ly;
       errmax(ti, 2,i) = max(max(abs(eomo1(i,:,:))));
     end
  end % end of t loop
  %plot(tOut, errmax(:,1,1)./errmax(:,2,1), tOut, errmax(:,1,2)./errmax(:,2,2), tOut, errmax(:,1,3)./errmax(:,2,3))
  plot(tOut, errmax(:,2,1), tOut, errmax(:,2,2), tOut, errmax(:,2,3))
  %plot(tOut, errmax(:,1,1), tOut, errmax(:,1,2), tOut, errmax(:,1,3))
  end

  tic
  T0UU = zeros(size(uOut, 1), params.Nx, params.Ny, 3, 3); 
  T1UU = zeros(size(uOut, 1), params.Nx, params.Ny, 3, 3); 
  T2UU = zeros(size(uOut, 1), params.Nx, params.Ny, 3, 3); 
  for ti = 1:size(uOut,1)
    [T0UU(ti, :,:,:,:) T1UU(ti,:,:,:,:) T2UU(ti,:,:,:,:)] = h2stressenergy(params, uOut(ti, :), udtOut(ti, :), udtdtOut(ti,:));
  end
  disp(sprintf('computed second-order stress energy tensor in %.2f sec', toc));
  size(uOut,1)

  tlen = size(uOut,1)-derorder;
  T0UUd = zeros(tlen, Nx, Ny, 3, 3, 3);
  T1UUd = zeros(tlen, Nx, Ny, 3, 3, 3);
  T2UUd = zeros(tlen, Nx, Ny, 3, 3, 3);
  T0UUd(:,:,:,:,:,1) = reshape(conv2(reshape(T0UU,size(uOut,1),Nx.*Ny.*3.*3), derker, 'valid')./params.Dt, tlen, Nx, Ny, 3, 3, 1);
  T1UUd(:,:,:,:,:,1) = reshape(conv2(reshape(T1UU,size(uOut,1),Nx.*Ny.*3.*3), derker, 'valid')./params.Dt, tlen, Nx, Ny, 3, 3, 1);
  T2UUd(:,:,:,:,:,1) = reshape(conv2(reshape(T2UU,size(uOut,1),Nx.*Ny.*3.*3), derker, 'valid')./params.Dt, tlen, Nx, Ny, 3, 3, 1);

  for ti=1:size(uOut,1)
    for i=1:3
      for j=1:3
        T0UUd(ti,:,:,i,j,2) = dx(reshape(T0UU(ti,:,:,i,j), Nx,Ny));
        T0UUd(ti,:,:,i,j,3) = dy(reshape(T0UU(ti,:,:,i,j), Nx,Ny));
        T1UUd(ti,:,:,i,j,2) = dx(reshape(T1UU(ti,:,:,i,j), Nx,Ny));
        T1UUd(ti,:,:,i,j,3) = dy(reshape(T1UU(ti,:,:,i,j), Nx,Ny));
        T2UUd(ti,:,:,i,j,2) = dx(reshape(T2UU(ti,:,:,i,j), Nx,Ny));
        T2UUd(ti,:,:,i,j,3) = dy(reshape(T2UU(ti,:,:,i,j), Nx,Ny));
      end
    end
  end

  divT0U = zeros(tlen, Nx, Ny, 3);
  divT1U = zeros(tlen, Nx, Ny, 3);
  divT2U = zeros(tlen, Nx, Ny, 3);

  for ti=1:tlen
    for xi=1:Nx
      for yi=1:Ny
        for i=1:3
          divT0U(ti,xi,yi,:) = divT0U(ti,xi,yi,:) + reshape(T0UUd(ti,xi,yi,:,i,i), 1, 1, 1, 3);
          divT1U(ti,xi,yi,:) = divT1U(ti,xi,yi,:) + reshape(T1UUd(ti,xi,yi,:,i,i), 1, 1, 1, 3);
          divT2U(ti,xi,yi,:) = divT2U(ti,xi,yi,:) + reshape(T2UUd(ti,xi,yi,:,i,i), 1, 1, 1, 3);
        end
      end
    end
  end

  divTUint = zeros(tlen, 3);
  divT0Uint = zeros(tlen, 3);
  divT1Uint = zeros(tlen, 3);
  divT2Uint = zeros(tlen, 3);

  divTUint = reshape(sum(sum(abs(divT0U+divT1U+divT2U),2),3),tlen,3);
  divT0Uint = reshape(sum(sum(abs(divT0U),2),3),tlen,3);
  divT1Uint = reshape(sum(sum(abs(divT1U),2),3),tlen,3);
  divT2Uint = reshape(sum(sum(abs(divT2U),2),3),tlen,3);

  deli = [1:(derorder/2) (length(tOut)-derorder/2+1):length(tOut)];
  tpOut = tOut;
  tpOut(deli, :) = [];

hold off
  subplot(1,1,1)
  plot(tpOut, divTUint(:,1), 'r', tpOut, divTUint(:,2), 'r', tpOut, divTUint(:,3), 'r');
  hold on
  %subplot(1,4,2)
  plot(tpOut, divT0Uint(:,1), 'b', tpOut, divT0Uint(:,2), 'b', tpOut, divT0Uint(:,3), 'b');
  %subplot(1,4,3)
  plot(tpOut, divT1Uint(:,1), 'g', tpOut, divT1Uint(:,2), 'g', tpOut, divT1Uint(:,3), 'g');
  %subplot(1,4,4)
  plot(tpOut, divT2Uint(:,1), 'm', tpOut, divT2Uint(:,2), 'm', tpOut, divT2Uint(:,3), 'm');




end
